Introduction

Soybean (Glycine max [L.] Merr.) is a primary sources of plant protein and edible oil worldwide, with seeds rich in protein (about 40%) and oil (about 20%) (Chiari et al. 2004). Soybean seed protein content (PC) and oil content (OC) is quantitative traits influenced by both genetic and environmental factors (Liang et al. 2010). The genetic effects include additive effects, epistasis and interactions of quantitative trait loci (QTLs) with the environment. In the wake of the improvements in molecular technology and statistical methods in recent decades, many QTLs have been identified in crop species. Numerous studies have identified QTLs for pairs of traits (Brummer et al. 1997; Orf et al. 1999; Csanádi et al. 2001; Liang et al. 2010; Pathan et al. 2013; Wang et al. 2014a; Warrington et al. 2015; Qi et al. 2017); however, only a few such QTLs have been identified in multiple environments and multiple genetic backgrounds. For example, Brummer et al. (1997) identified QTLs for PC and OC in 8 soybean populations that were sensitive to environmental and genetic background; fewer than 15 stable QTLs were identified for each trait, and no population had more than 3 stable QTLs. Moreover, for OC, no stable QTLs were identified in 2 of the 8, although the other 6 populations each contained at least a single stable QTL, and one population had 3; for PC, at least one stable QTL was found in 8 populations. Orf et al. (1999) used amplified fragment length polymorphism (AFLP) markers and simple sequence repeat (SSR) markers in the three RIL populations derived from 3 parents, Minsoy, Noir 1 and Archer, in four environments. Five PC and 6 OC QTLs were detected, but most were identified in only one population, and no identical QTLs were identified in multiple populations. Wang et al. (2014a) detected 3-trait QTLs using 2 RIL populations in multiple environments, among which 9 PC and 8 OC QTLs were further confirmed by comparison with previously reported QTLs, and the other 8 were newly identified. Using MAS, a trait can be successfully expressed in a plant if the control of the related QTL is not affected by the environment or the genetic background.

In addition to additive effects, epistasis (additive × additive interaction) is another major genetic basis for complex phenotypic traits, playing a vital role in heterosis, breeding inhibition, adaptability, reproductive isolation and speciation (Yang and Zhu 2005). Many additive × additive (AA) and interaction with environment (AAE) epistatic QTLs for soybean have been detected in recent years (Hou et al. 2014; Qi et al. 2014; Wang et al. 2015a; Qi et al. 2017; Teng et al. 2017; Tan et al. 2018). Hou et al. (2014) mapped PC and OC QTLs using SSR markers derived from the strains Charleston and Dongnong594 and detected 3 epistatic-effect QTL pairs related to PC and 4 for OC; Qi et al. (2014) identified additive- and epistatic-effect QTLs for PC and OC in multiple environments in the same populations. Teng et al. (2017) detected 7 additive QTL pairs and 5 epistatic-effect QTL pairs for soybean seed oil quality. In summary, the identification of epistatic QTL interactions has largely been conducted using only separate single populations and separate environments, without consideration for the stability of the associations in multiple genetic backgrounds or environments.

In this study, we used two soybean RIL populations derived from the crosses Dongnong L13 × Henong 60 and Dongnong L13 × Heihe 36 and planted in 8 environments to identify AA and AAE QTLs for soybean seed PC and OC by ICIM and MCIM, with the goals of exploring the genetic architecture of PC and OC and improving the efficiency of MAS for soybean quality traits.

Materials and Methods

Plant materials and field design

Two populations, RIL3613 (Dongnong L13 × Heihe 36) and RIL6013 (Dongnong L13 × Henong 60), containing 134 and 156 RIL_2:8, respectively, were obtained from crosses between three soybean parents with major differences in quality and other characteristics, Dongnong L13 (PC 45.50%, OC 18.74%), Henong 60 (PC 38.47%, OC 22.25%), and Heihe 36 (PC 39.80%, OC 19.28%). Starting in the F₂ generation, the seeds of each single plant were propagated by single-seed descent, and RIL populations obtained after five successive generations of self-crossing in 2008 in Harbin (HRB; 45°75ʹ N, 126°63ʹ E), Heilongjiang, China, and Yacheng (17°50ʹ N, 109°00ʹ E) in Hainan Province, China, were used for map construction.

The parental lines and RILs were planted in 8 environments: in Keshan (KS; 48°25ʹ N, 125°64ʹ E) in 2013; in Harbin (HRB; 45°75ʹ N, 126°63ʹ E) in 2014; in Harbin and Keshan in 2015; in Acheng (AC; 45°52ʹ N, 126°95ʹ E), Shuangcheng (SC; 45°53ʹ N, 126°32ʹ E) and Harbin in 2016; and in Shuangcheng (SC; 45°53ʹ N, 126°32ʹ E) in 2017. Three replicate plantings of each line were grown in a randomized complete block design, using rows 3 m in length, 0.70 m apart, with the seeds in each individual row sown at 0.06-m intervals.

Measurement of oil and protein contents

Seed phenotypic measurements were obtained from ten mature plants randomly selected in the middle row of each plot. The PC and OCs of seed were determined three times with an Infratec 1241 Grain Analyzer (FOSS, Sweden) at the 13% moisture basis.

Variation analysis and heritability of phenotypic data

The significance of the differences in PC and OC between the two parents of each population was determined by Student's t test, and the significance of the genotype differences between RILs and environments was determined by ANOVA. The frequency distributions were analyzed with Microsoft Excel 2007. The following formulas were used to estimate heritability.

For single environments:

For the multi-environment average values:

Where h² is broad-sense heritability, is the variance of genotype, is the variance of error, indicates variance of genotype by environment effect, r is the number of replications and e is the number of environments in the study. , and were estimated using a mixed method implemented by Proc Mixed in SAS9.1 (SAS Institute Inc., USA).

QTL mapping

On account of the SSR linkage map constructed in the previous study (Ning et al. 2018). The total SSR linkage map lengths were 2849.54 cM and 1886.8 cM and the mean interval lengths were 21.92 cM and 16.13 cM for RIL3613 and RIL6013, respectively. The average of the quality traits for each strain was analyzed conjointly in multiple environments by the inclusive composite interval-mapping (ICIM) method (Li et al. 2006) and by composite interval mapping based on mixed linear models (MCIM) (Yang et al. 2008). Using the software QTL IciMapping v4.2, the ICIM-ADD and ICIM-EPI algorithms of the MET model of ICIM were applied to analyze the additive-effect and epistatic-effect QTLs. The mapping step was set to 2.0 cM, and LOD thresholds were determined by 1000 permutation tests combining probability of 0.05 for type I error. QTL Network 2.0 software was used to detect additive- and epistatic-effect QTLs based on MCIM. One- and two-dimensional genome scans for QTLs were performed using a 10-cM testing window, a 0.1 cM walk speed and a 0.5 cM filtration window size. To control the experimental type I error rate, a critical F value using the Satterthwaite method was estimated by performing a permutation test 1,000 times. The naming of QTLs followed the QTL nomenclature described by McCouch et al. (1997).

Results

Phenotypic variation

To investigate the genetic basis for soybean seed protein content (PC) and oil content (OC); we assessed PC and OC in soybeans from two RIL populations in eight different environments (defined year and location). The data revealed significant variation among both the RIL3613 and the RIL6013 lines (Table 1, 2 and 3); the minimum and maximum values differed widely, the skewness and kurtosis values were <1.00 and the data were normally distributed (Fig. 1). An ANOVA to detect the interactions of PC and OC with genotype, with environment and with genotype × environment showed significant interactions (P < 0.05).

Additive effect QTLs

In this study, we identified a total of 33 and 41 QTLs related to the two traits, located on 18 of the 20 soybean chromosomes (all but K and N), in the RIL3613 and RIL6013 populations, respectively, grown under the eight environments (Fig. 2).

Fig. 1: Frequency distribution of protein and oil contents in two populations under eight environments

E1: Keshan in 2013; E2: Harbin in 2014; E3: Harbin in 2015; E4: Keshan in 2015; E5: Acheng in 2016; E6: Shuangcheng in 2016; E7: Harbin in 2016; E8: Shuangcheng in 2017

In the RIL3613 population, we identified 30 PC and 3 OC additive-effect QTLs in the 17 soybean linkage group (barring K, L and N); the LOD values ranged from 2.53 to 7.88 and from 6.79 to 11.23 for PC and OC QTLs, respectively, and the proportion of phenotypic variability explained (PVE) values were 2.54–13.88% and 13.8–38.44%, respectively (Fig. 2 and Table 4). qPro-D2-3, qOil-A2-1 and qOil-G-1 had PVE values of more than 10%. Ten of the QTLs for PC (qPro-A2-1, qPro-B1-1, qPro-C1-3, qPro-D1a-3, qPro-G-6, qPro-H-1, qPro-I-1, qPro-J-3, qPro-L-2 and qPro-O-1) had positive additive effects, meaning that the alleles derived from Dongnong L13 increased PC (ADD > 0.1). Nine PC QTLs (qPro-A2-2, qPro-D1a-2, qPro-D1b-1, qPro-D1b-5, qPro-D2-3, qPro-F-4, qPro-G-1, qPro-J-2 and qPro-L-4) and one 1 OC QTL (qOil-G-3) had negative additive effects, with the alleles from Heihe 36 increasing PC or OC (ADD < –0.1%).

In the RIL6013 population, we identified 21 PC and 20 OC additive-effect QTLs on 16 soybean chromosomes (linkage groups A1, B1, B2, C1, C2, D1a, D1, D2, E, F, G, H, I, J, M and O); the LOD values ranged from 2.53 to 4.83 and 2.52 to 6.53, respectively, with PVEs of 2.74–11.64% and 2.99–7.96% (Fig. 2 and Table 5). Moreover, the PVEs of qPro-E-1, qPro-F-6, qPro-M-5, qOil-C1-1, qOil-D1a-2, qOil-D2-1, qOil-H-1 and qOil-I-2 were all more than 10%. Five QTLs for PC (qPro-A1-3, qPro-C1-2, qPro-D1a-1, qPro-G-4 and qPro-G-5) and one QTL for OC (qOil-D1b-1) had positive additive effects, meaning that the alleles derived from Dongnong L13 enhanced the PC or OC (ADD > 0.1%), while 2 PC QTLs (qPro-E-1 and qPro-F-6) and two OC QTLs (qOil-D1a-2 and qOil-D2-2) had negative additive effects, with the alleles from Henong 60 increasing the PC or OC (ADD < –0.1).

Seven QTLs were detected by both methods (Table 4 and 5); among these, qPro-G-3, qPro-G-6 and qPro-C1-1 had positive additive effects, meaning that the alleles from Dongnong L13 enhanced PC, whereas qPro-D2-3, qOil-A2-1 and qOil-H-1 had negative additive effects, with the alleles from Dongnong L13 reducing OC.

A total of seven QTLs with multiple effects simultaneously controlled PC and OC. Among these, the QTL qPro-D1b-3 (Satt041-Satt546, 84.04–87.19 cM) for PC was found in both the RIL3613 and RIL6013 populations, and had a positive additive effect, indicating that the allele from Dongnong L13 increased PC. Meanwhile, six QTL SSR intervals (Satt276-Sat_171, Sct_067-Satt589, Sat_413-Sat_160, Satt685-Satt231, AZ254740-Satt570, Satt414-Sat_255) simultaneously control PC and OC with opposite additive effects, which implies that it may be difficult to improve PC and OC at the same time through the use of these QTLs.

Epistatic-effect QTLs

We identified 18 epistatic-effect QTL pairs for either PC or OC in the two RIL populations under eight environments by multiple-environment interaction (AAE) analysis using ICIM and MCIM methods for the combinations (Table 6 and Fig. 3). Among them, nine epistatic-effect QTL pairs related to PC and two pairs related to OC had positive additive effects and the other four pairs for PC and three pairs for OC had negative additive effects.

Table 1: Summarization of protein content in eight environments

Environment ^A	Parents		RILs						F	h^2^B
	Dongnong L13	Heihe 36	Average	Std	Min	Max	Kurtosis	Skewness
RIL3613
2013KS	41.81	41.00	43.39	2.07	38.00	47.37	-0.56	-0.34	229.41^**C	0.987
2014HRB	40.60	40.00	42.35	2.23	36.74	46.14	-0.61	-0.35	245.90^**	0.988
2015HRB	43.83	41.50	42.62	1.26	38.60	44.80	-0.30	-0.51	82.08^**	0.964
2015KS	44.10	41.80	41.71	1.44	37.59	44.77	-0.26	-0.22	113.12^**	0.974
2016AC	44.20	41.90	41.35	1.27	37.80	44.20	-0.17	-0.27	81.90^**	0.964
2016SC	43.20	41.70	41.99	1.31	37.80	45.00	0.79	-0.80	102.85^**	0.971
2016HRB	43.40	41.00	41.51	1.24	37.40	44.00	0.40	-0.71	86.24^**	0.966
2017SC	40.50	44.00	42.50	0.99	38.10	44.30	1.93	-0.95	57.23^**	0.949
RIL6013	Dongnong L13	Henong 60
2013KS	40.70	43.50	44.13	1.90	39.18	48.30	-0.44	-0.03	359.91^**	0.992
2014HRB	41.20	42.30	43.63	1.55	39.63	47.49	-0.41	0.21	232.82^**	0.987
2015HRB	41.50	42.20	43.58	0.92	39.90	46.00	1.84	-0.71	81.79^**	0.964
2015KS	40.90	43.10	42.78	1.14	39.21	45.30	0.38	-0.30	156.12^**	0.981
2016AC	40.80	42.20	42.36	1.10	38.90	45.70	1.03	0.18	128.51^**	0.977
2016SC	41.20	43.60	42.93	1.12	39.40	45.80	0.50	-0.57	141.90^**	0.979
2016HRB	41.70	42.20	42.18	1.06	39.30	45.30	0.19	-0.32	124.94^**	0.976
2017SC	43.50	43.00	42.77	0.99	39.50	44.60	0.42	-0.70	111.46^**	0.974

A: 2013KS means Keshan in 2013; 2014HRB means Harbin in 2014; 2015 HRB means Harbin in 2015; 2015KS means Keshan in 2015; 2016AC means Acheng in 2016; 2016SC means Shuangcheng in 2016; 2016HRB means Harbin in 2016. 2017SC means Shuangcheng in 2017

B: h² means broad-sense heritability

C: ** means significant at 0.01 levels

Table 2: Summarization of oil content in eight environments

Environment ^A	Parents		RILs						F	h^2B
	Dongnong L13	Heihe 36	Average	Std	Min	Max	Kurtosis	Skewness
RIL3613
2013KS	19.90	20.10	17.07	0.86	15.19	18.89	-0.69	0.15	19.98^**C	0.864
2014HRB	20.10	20.30	20.51	1.01	17.89	22.37	-0.29	-0.36	27.65^**	0.899
2015HRB	19.28	19.80	20.43	0.51	18.54	21.92	0.93	-0.17	6.43^**	0.644
2015KS	20.25	20.15	19.69	0.85	16.72	21.83	0.23	-0.10	17.82^**	0.849
2016AC	19.28	18.88	20.26	0.59	18.81	22.25	0.73	-0.27	10.13^**	0.753
2016SC	19.98	20.53	20.18	0.61	18.51	21.78	0.29	-0.28	11.53^**	0.778
2016HRB	20.26	20.90	20.68	0.58	18.39	22.10	2.53	-0.79	10.58^**	0.761
2017SC	21.90	21.50	21.28	0.40	20.10	22.40	0.79	-0.25	6.16^**	0.632
RIL6013	Dongnong L13	Henong 60
2013KS	20.71	20.18	17.03	0.90	14.40	20.13	0.50	0.39	50.04^**	0.942
2014HRB	20.87	20.22	19.52	0.93	16.73	22.39	0.47	0.10	48.28^**	0.940
2015HRB	21.21	20.92	20.16	0.50	18.45	22.63	4.86	0.56	13.63^**	0.808
2015KS	20.07	20.78	19.34	0.82	17.19	21.88	0.06	-0.14	35.69^**	0.920
2016AC	20.09	20.38	20.13	0.54	18.16	21.14	0.83	-0.87	16.55^**	0.838
2016SC	20.94	20.24	20.00	0.54	18.39	21.21	0.20	-0.48	16.30^**	0.836
2016HRB	20.93	20.78	20.68	0.40	18.98	21.74	1.93	-0.81	9.93^**	0.748
2017SC	21.40	21.40	21.32	0.38	19.20	22.30	6.08	-1.34	8.41^**	0.712

A: 2013KS means Keshan in 2013; 2014HRB means Harbin in 2014; 2015 HRB means Harbin in 2015; 2015KS means Keshan in 2015; 2016AC means Acheng in 2016; 2016SC means Shuangcheng in 2016; 2016HRB means Harbin in 2016; 2017SC means Shuangcheng in 2017

B: h² means broad-sense heritabili

C: ** means significant at 0.01 levels

Table 3: Analysis of variance and heritability on protein and oil contents across multiple environments

Population	Trait^A	Max-imum	Min-imum	Mean	Standard deviation	CV	FE^B	FG^C	FG×E^D	h^2^E
RIL3613	PC	47.37	36.74	42.18	1.64	3.34	3196.89^**F	242.2^**	107.10^**	0.572
	OC	22.40	15.19	19.99	1.40	3.32	5370.70^**	19.13^**	12.45^**	0.369
RIL6013	PC	48.30	39.18	43.06	1.41	2.86	7062.22^**	227.58^**	161.06^**	0.311
	OC	22.63	15.17	14.40	1.39	3.26	13450.0^**	33.09^**	23.30^**	0.317

A: PC means protein content; OC means oil content

B: F _E means F value for environment effects

C: F _G means F value for genetic effects

D: F _{G × E} means F value for genotype × environment interaction effects

E: h² means broad-sense heritability

F: ** means significant at 0.01 levels

We detected 13 sites of pairwise interaction related to PC by AA and AAE analysis in the two RIL populations (Table 6 and Fig. 3). The AA values ranged from 0.79% to 2.72%, the PVEs for AA ranged from 0.07 to 3.09%, the total PVE for AAE was 12.87%, the PVEs for AAE interaction ranged from 0.79 to 2.72%, and the total PVE

Fig. 2: QTL associated with protein (red bars) and oil (blue bars) contents in RIL3613 and RIL6013
Table 4: Additive QTLs associated with protein and oil contents in RIL3613

QTL	Chr	Marker interval	Region in public map	Analysis method	LOD^A	PVE (%)^B	h^2(%)^C	ADD^D	AE1^E	AE2^E	AE3^E	AE4^E	AE5^E	AE6^E	AE7^E	AE8^E
PC
qPro-A1-2	A1	Satt717-Sat_171	51.95-57.79	ICIM	2.85	3.70		-0.09	-0.219	0.141	0.032	0.114	0.140	-0.017	-0.173	-0.018
qPro-A2-1	A2	Sct_067-Satt589	14.99-33.95	ICIM	4.16	7.82		0.10	0.129	0.422	0.045	-0.087	-0.308	-0.101	0.001	-0.101
qPro-A2-2	A2	Satt424-Satt233	60.59-100.08	ICIM	4.74	5.83		-0.17	-0.044	-0.116	0.074	-0.006	0.003	-0.098	-0.037	0.225
qPro-B1-1	B1	Satt197-Satt359	46.38-102.55	ICIM	4.11	5.29		0.12	0.237	-0.085	0.181	0.069	-0.199	-0.048	-0.111	-0.044
qPro-C1-3	C1	Sat_140-Sat_416	41.43-76.41	ICIM	5.36	7.18		0.16	0.242	-0.055	0.042	0.164	-0.032	-0.063	-0.250	-0.048
qPro-C2-1	C2	Sat_336-Satt681	3.15-51.84	ICIM	2.91	3.59		0.08	0.151	0.030	0.124	0.148	-0.194	-0.177	-0.033	-0.049
qPro-C2-5	C2	Satt202-Satt316	126.23-127.66	ICIM	2.53	4.20		-0.05	0.244	-0.300	0.091	0.061	-0.017	-0.218	0.058	0.082
qPro-D1a-2	D1a	Sat_346-Satt515	53.66-55.68	ICIM	4.08	4.34		-0.10	0.130	-0.216	0.258	0.225	-0.100	-0.160	0.022	-0.158
qPro-D1a-3	D1a	Sa_346-Satt198	53.66-68.62	ICIM	5.27	4.74		0.28	-0.299	-0.209	-0.230	-0.284	0.352	0.268	0.351	0.051
qPro-D1b-1	D1b	Satt698-AI856415	38.04-50.11	ICIM	3.21	3.49		-0.15	-0.070	0.372	-0.248	-0.124	0.049	-0.043	-0.108	0.172
qPro-D1b-3	D1b	Satt041-Satt546	84.04-87.19	ICIM	3.48	3.12		-0.03	0.003	-0.172	-0.001	0.104	-0.163	0.048	-0.114	0.295
qPro-D1b-5	D1b	Sat_069-Satt271	102.59-137.05	ICIM	3.21	5.16		-0.13	-0.045	-0.275	-0.020	-0.024	0.008	0.011	0.144	0.200
qPro-D2-3	D2	Sat_001-Sat_326	92.12-112.84	ICIM	7.23	13.88		-0.17	-0.510	-0.278	0.245	0.258	0.022	-0.003	0.040	0.226
				MCIM			1.41	-0.25	-0.369	-0.586	0.338	0.245	0.038	-0.014	0.049	0.304
qPro-F-1	F	GMRUBP-Sat_262	0-9.69	ICIM	2.73	2.54		0.05	-0.028	-0.248	0.167	-0.150	0.013	0.117	-0.019	0.147
qPro-F-4	F	Sat_039-SOYHSP176	27.87-68.44	ICIM	3.36	3.81		-0.10	0.015	-0.148	0.045	0.060	-0.073	-0.087	-0.077	0.265
qPro-G-1	G	Sat_210-Satt688	3.7-12.54	ICIM	3.39	5.34		-0.12	-0.172	-0.275	0.074	0.177	-0.113	0.025	0.038	0.246
qPro-G-3	G	AZ254740-Satt570	8.23-12.74	ICIM	7.88	8.26		0.00	-0.147	-0.321	0.404	0.316	-0.213	-0.075	0.134	-0.099
				MCIM			0.14	0.05	-0.118	-0.236	0.284	0.347	-0.157	-0.073	0.086	-0.127
qPro-G-6	G	Satt503-Satt288	68.76-76.76	ICIM	6.67	7.17		0.13	0.102	0.007	0.115	0.336	-0.205	-0.143	-0.054	-0.158
				MCIM			0.99	0.35	0.477	-0.092	0.428	0.655	-0.400	-0.468	-0.123	-0.465
qPro-H-1	H	Sat_200-Satt353	3.02-8.48	ICIM	3.24	3.71		0.14	0.018	-0.665	0.158	0.158	-0.010	0.307	-0.204	0.239
qPro-I-1	I	Satt367-Satt270	27.98-50.11	ICIM	4.55	8.78		0.16	0.453	0.090	-0.065	-0.074	-0.201	-0.092	-0.056	-0.055
qPro-I-2	I	Satt354-Sct_189	46.22-113.76	ICIM	2.74	2.89		0.06	0.005	0.197	0.178	0.171	-0.213	-0.177	-0.004	-0.158
qPro-J-2	J	Satt414-Sat_350	37.04-55.73	ICIM	6.11	8.84		-0.21	-0.285	-0.145	-0.021	-0.062	-0.047	0.186	0.063	0.311
qPro-J-3	J	Satt654-Sat_224	38.09-75.12	ICIM	3.27	4.83		0.18	0.181	0.135	-0.041	0.139	-0.223	0.041	-0.167	-0.066
qPro-J-4	J	Sct_193-Satt183	41.5-42.5	ICIM	3.62	5.11		0.05	0.360	0.051	0.094	-0.050	-0.072	-0.157	0.032	-0.259
qPro-L-1	L	Satt182-Sat_134	14.03-28.27	ICIM	3.15	5.38		-0.07	0.067	-0.381	0.089	0.198	-0.125	-0.067	0.087	0.132
qPro-L-2	L	Sat_134-Sat_191	28.27-32	ICIM	6.61	9.38		0.17	0.305	0.202	0.072	0.031	-0.251	-0.133	-0.108	-0.119
				MCIM			1.99	0.45	0.467	0.398	0.055	0.038	-0.270	-0.192	-0.155	-0.334
qPro-L-4	L	Sat_099-Satt229	78.23-93.88	ICIM	3.92	4.50		-0.22	-0.323	0.088	0.231	0.159	-0.390	0.126	0.056	0.054
qPro-M-3	M	Satt567-Satt697	33.47-85.34	ICIM	5.01	4.78		0.00	-0.123	-0.038	-0.207	-0.218	0.217	0.297	-0.049	0.121
qPro-M-4	M	Sat_121-Satt346	103.98-112.79	ICIM	5.99	6.85		0.08	-0.176	0.104	0.189	0.380	-0.134	-0.168	-0.015	-0.182
qPro-O-1	O	Satt358-Sat_303	5.44-20.93	ICIM	4.22	8.76		0.18	0.040	0.598	-0.201	-0.077	-0.202	0.047	-0.153	-0.052
OC
qOil-A2-1	A2	Sct_067-Satt589	14.99-33.95	ICIM	11.23	38.44		-0.08	-0.262	-0.260	0.022	0.050	0.114	0.074	0.127	0.134
				MCIM			0.42	-0.10	-0.238	-0.271	0.028	0.026	0.119	0.077	0.127	0.133
qOil-G-1	G	AZ254740-Satt570	8.23-12.74	ICIM	6.79	13.80		-0.04	0.146	0.016	-0.107	-0.204	0.100	0.004	0.044	0.001
qOil-G-3	G	Satt503-Satt288	68.76-76.76	MCIM			0.45	-0.15	-0.142	-0.054	-0.005	-0.051	0.083	0.038	0.024	0.106

A: LOD, log of odd

B: PVE means phenotypic variation explanation ration

C: h² means phenotypic variation explained by additive QTL

D: ADD means additive effects

E: Additive by environment interaction effect. E1: Keshan in 2013; E2: Harbin in 2014; E3: Harbin in 2015; E4: Keshan in 2015; E5: Acheng in 2016; E6: Shuangcheng in 2016; E7: Harbin in 2016; E8: Shuangcheng in 2017

Table 5: Additive QTLs associated with protein and oil contents in RIL6013

QTL	Chr	Marker interval	Region in public map	Analysis method	LOD^A	PVE (%)^B	h^2(%)^C	ADD^D	AE1^E	AE2^E	AE3^E	AE4^E	AE5^E	AE6^E	AE7^E	AE8^E
PC
qPro-A1-1	A1	Satt276-Sat_171	17.16-57.79	ICIM	2.70	5.92		0.08	0.190	0.028	0.033	0.065	-0.143	-0.086	0.015	-0.102
qPro-A1-3	A1	Satt545-Satt200	71.38-92.88	ICIM	3.68	4.04		0.14	-0.090	-0.227	0.168	0.196	0.111	0.074	-0.021	-0.211
qPro-B1-2	B1	Sat_128-Sat_095	53.41-81.3	ICIM	2.53	2.74		-0.06	0.122	0.026	-0.097	-0.114	-0.008	-0.003	0.037	0.039
qPro-B2-1	B2	Sat_230-Satt474	72.08-75.34	ICIM	3.19	4.61		-0.08	0.193	0.145	-0.047	-0.105	-0.122	-0.058	-0.039	0.034
qPro-C1-1	C1	Satt565-Satt713	0-88.94	ICIM	3.04	9.04		0.05	0.360	0.089	-0.137	-0.191	-0.010	0.021	-0.050	-0.083
				MCIM			0.35	0.11	0.384	0.286	-0.138	-0.213	-0.065	-0.033	-0.114	-0.107
qPro-C1-2	C1	Sat_367-Sat_140	28.04-41.43	ICIM	2.73	3.86		0.11	-0.001	0.065	0.102	0.100	-0.005	-0.064	-0.107	-0.091
qPro-C2-3	C2	Satt376-Satt307	97.83-121.26	ICIM	4.24	7.88		-0.01	0.466	-0.153	0.131	0.166	-0.118	-0.134	-0.361	0.002
qPro-C2-4	C2	Satt277-Satt316	107.58-127.66	ICIM	2.9	3.83		-0.04	0.063	-0.038	-0.098	-0.108	-0.106	0.015	0.148	0.123
qPro-D1a-1	D1a	Sat_413-Sat_160	5.93-104.27	ICIM	2.56	4.98		0.17	0.092	0.328	-0.160	-0.268	0.079	0.018	-0.014	-0.074
qPro-D1b-3	D1b	Satt041-Satt546	84.04-87.19	ICIM	2.76	5.11		0.02	0.180	0.160	-0.163	-0.123	-0.022	-0.140	0.105	0.002
qPro-D2-1	D2	Satt154-Satt669	57.07-67.7	ICIM	3.46	8.69		0.06	0.394	-0.072	-0.067	-0.098	-0.138	0.042	-0.189	0.127
qPro-D2-2	D2	Sat_194-Sat_001	86.69-92.12	ICIM	2.55	4.73		0.02	0.224	0.090	-0.142	-0.141	-0.040	0.078	0.073	-0.142
qPro-E-1	E	Satt483-Satt553	44.98-67.91	ICIM	2.97	10.02		-0.10	-0.332	-0.061	0.016	0.041	0.043	0.124	0.045	0.123
qPro-E-2	E	Satt685-Satt231	56.7-70.23	ICIM	3.05	6.36		0.01	0.305	-0.006	0.011	0.014	-0.096	-0.103	0.109	-0.233
qPro-F-6	F	Satt334-Sat_417	78.05-135.94	ICIM	4.83	11.64		-0.10	-0.130	-0.334	0.190	0.168	0.018	0.043	-0.041	0.087
qPro-G-4	G	Satt570-AW734137	12.74-15.63	ICIM	3.16	6.16		0.11	0.078	0.143	0.030	0.097	-0.151	-0.093	-0.004	-0.099
qPro-G-5	G	Satt352-Satt564	50.52-57.32	ICIM	3.89	5.88		0.27	-0.183	0.366	0.023	-0.102	-0.209	-0.135	0.605	-0.366
qPro-I-3	I	Sat_268-Sat_170	55.09-75	ICIM	2.95	9.34		0.00	0.454	-0.196	0.003	-0.063	-0.055	0.050	-0.035	-0.158
qPro-J-1	J	Satt414-Sat_255	37.04-43.84	ICIM	3.92	8.28		-0.08	-0.195	-0.160	0.195	0.176	-0.027	-0.072	-0.054	0.137
qPro-M-5	M	Satt210-Satt346	112.08-112.79	ICIM	4.52	11.22		0.06	0.429	-0.102	-0.181	-0.208	-0.068	0.085	-0.024	0.068
qPro-O-2	O	Sat_303-Satt633	20.93-56.93	ICIM	2.61	8.78		0.07	0.376	0.056	-0.146	-0.091	-0.035	-0.088	0.003	-0.074

Table 5: Continued

for environmental interaction was 19.77% for PC. Four pairwise interaction sites, qPro-D1b-4~qPro-N-1, qPro-M-2~qPro-F-2, qPro-H-3~qPro-G-2 and qPro-B2~2-qPro-J-7, showed negative epistatic effects, while the Table 5: Continued

OC
qOil-A1-1	A1	Satt276-Sat_171	17.16-57.79	ICIM	3.00	4.09		-0.02	0.085	-0.022	0.037	0.030	-0.074	-0.051	-0.020	0.015
qOil-A1-2	A1	Satt545-Satt174	71.38-88.58	ICIM	2.69	2.99		-0.03	0.143	0.121	0.001	-0.006	0.019	-0.075	-0.099	-0.105
qOil-A1-3	A1	Sat_267-Satt200	78.44-92.88	ICIM	3.04	4.19		-0.05	0.035	0.029	-0.058	-0.005	0.056	-0.086	0.002	0.027
					2.52	4.40		-0.00	0.095	-0.057	0.056	0.064	-0.089	-0.056	0.000	-0.013
qOil-B1-1	B1	Satt197-Sat_123	46.38-100.88	ICIM	2.70	4.97		-0.02	-0.110	0.070	-0.032	0.094	-0.031	-0.058	-0.001	0.068
					5.63	8.23		0.04	-0.038	-0.421	-0.068	-0.087	0.033	0.278	0.160	0.142
qOil-B2-1	B2	Satt168-Sat_009	55.2-78.66	ICIM	2.72	4.59		0.02	-0.021	-0.163	0.027	0.086	0.003	0.006	0.027	0.036
qOil-C1-1	C1	Satt396-Sat_367	24.11-28.04	ICIM	4.24	10.68		-0.04	-0.239	0.057	-0.022	-0.013	0.088	0.000	0.029	0.100
qOil-C2-2	C2	Sat_246-Satt277	91.8-107.58	ICIM	2.88	3.61		0.00	-0.044	0.015	0.081	0.115	-0.038	-0.056	-0.082	0.009
qOil-D1a-2	D1a	Sat_413-Sat_160	5.93-104.27	ICIM	3.97	10.36		-0.11	0.138	-0.322	0.086	0.042	-0.020	0.021	0.017	0.038
qOil-D1b-1	D1b	Staga002-Sat_289	126.44-131.91	ICIM	4.48	7.35		0.10	-0.027	-0.033	0.120	0.129	0.009	0.058	-0.102	-0.154
qOil-D2-1	D2	Sat_333-Sat_194	5.83-86.69	ICIM	6.53	17.96		-0.07	-0.268	-0.065	0.112	0.064	0.022	0.012	0.044	0.080
qOil-D2-2	D2	Sat_194-Sat_001	86.69-92.12	MCIM			0.28	-0.11	-0.262	-0.062	0.097	0.025	0.028	0.040	0.060	0.074
qOil-E-1	E	Satt685-Satt231	56.7-70.23	ICIM	3.74	6.35		0.00	0.073	-0.100	0.003	-0.153	0.125	0.165	-0.035	-0.078
qOil-F-1	F	Satt030-Sat_240	3.95-25.58	ICIM	3.73	8.57		-0.01	0.103	0.198	-0.035	-0.302	0.067	0.101	-0.010	-0.122
qOil-G-2	G	Satt688-Satt570	12.54-12.74	ICIM	2.61	5.99		-0.06	-0.034	-0.035	-0.008	-0.050	0.039	0.024	0.040	0.023
qOil-H-1	H	Satt181-Satt434	91.12-105.73	ICIM	6.25	15.7		-0.02	-0.211	-0.130	0.042	0.184	-0.031	0.011	0.057	0.078
				MCIM			0.07	-0.06	-0.252	-0.193	0.040	0.197	-0.010	0.037	0.096	0.084
qOil-I-1	I	Satt571-Satt367	18.5-27.98	ICIM	2.99	7.68		-0.06	0.025	0.000	-0.031	-0.157	0.055	0.012	0.054	0.043
qOil-I-2	I	Sat_170-Satt330	75-77.83	ICIM	6.08	12.28		0.03	0.175	0.036	0.016	0.125	-0.109	-0.089	-0.052	-0.101
qOil-J-1	J	Satt414-Sat_255	37.04-43.84	ICIM	3.56	6.81		0.05	0.045	0.066	-0.049	-0.090	0.050	0.025	0.001	-0.049
qOil-O-1	O	BF008905-Sat_221	28.95-51	ICIM	3.13	7.99		-0.05	0.200	0.264	-0.078	-0.294	-0.018	-0.036	-0.023	-0.015
qOil-O-2	O	Sat_221-Sat_341	51-67.93	ICIM	3.4	8.89		0.06	0.018	0.078	0.009	0.104	-0.009	-0.109	-0.074	-0.018

A: LOD, log of odd

B: PVE means phenotypic variation explanation ration

C: h² means phenotypic variation explained by additive QTL

D: ADD means additive effects

E: Additive by environment interaction effect. E1: Keshan in 2013; E2: Harbin in 2014; E3: Harbin in 2015; E4: Keshan in 2015; E5: Acheng in 2016; E6: Shuangcheng in 2016; E7: Harbin in 2016; E8: Shuangcheng in 2017

Table 6: Epistatic QTL for protein and oil contents

Trait	Popu-lation	QTL_i	Marker Interval	QTL_j	Marker Interval	Analysis method	AA^A	h^2(AA)^B (%)	h^2(AAE)^C (%)	AAE1^D	AAE2^D	AAE3^D	AAE4^D	AAE5^D	AAE6^D	AAE7^D	AAE8^D
PC	RIL3613	qPro-D1b-4	Sat_069Sat_183	qPro-N-1	Satt631-Satt125	ICIM	-0.115	0.53	1.84	-0.275	0.132	-0.279	-0.243	0.284	0.17	0.107	0.104
		qPro-D1a-4	Satt515-Satt254	qPro-M-4	Sat_121-Satt346	ICIM	0.191	1.11	1.51	0.191	0.219	0.281	0.166	-0.25	-0.342	-0.123	-0.142
		qPro-C2-2	Satt640-Satt281	qPro-M-4	Sat_121-Satt346	ICIM	0.167	1.02	2.14	0.496	0.076	0.129	0.043	-0.16	-0.311	-0.163	-0.109
		qPro-M-2	Sat_389Satt697	qPro-F-2	Satt030-Sat_262	ICIM	-0.227	1.81	0.79	-0.238	-0.054	-0.041	-0.126	0.101	0.09	0.034	0.234
		qPro-D1b-3	Satt041-Satt546	qPro-F-5	Satt510-Satt334	ICIM	0.224	1.77	1.3	0.371	0.154	-0.011	-0.028	-0.147	-0.367	0.027	-0.346
		qPro-J-5	Sct_193Sat_255	qPro-J-6	Sat_255-Satt620	ICIM	0.056	0.12	2.72	0.329	0.004	0.285	0.286	-0.222	-0.454	-0.164	-0.062
		qPro-H-3	Satt293-Satt434	qPro-G-2	Sat_210AW734137	ICIM	-0.223	1.79	1.39	-0.114	-0.295	-0.121	-0.208	0.15	0.158	0.098	0.332
		qPro-D1a-4	Satt515-Satt254	qPro-M-1	Sat_389-Satt245	MCIM	0.200	0.88	1.57	0.260	0.179	-0.035	-0.026	-0.212	-0.055	-0.178	0.064
		qPro-D1a-2	Sat_346Satt515	qPro-M-3	Satt567-Satt697	MCIM	0.543	3.09	1.15	0.112	0.035	-0.059	-0.075	0.044	0.035	-0.05	-0.042
		qPro-D1b-2	Satt698-Satt271	qPro-L-3	Satt497-Sat_099	MCIM	0.063	0.07	2.37	-0.326	0.663	-0.043	0.067	-0.217	-0.058	-0.223	0.145
		qPro-J-4	Sct_193Satt183	qPro-J-6	Sat_255-Satt620	MCIM	0.032	0.12	1.31	0.695	-0.162	0.445	0.458	-0.429	-0.794	-0.151	-0.062
	RIL6013	qPro-H-2	Satt293-Satt181	qPro-F-3	Satt030-Sat_240	MCIM	0.062	0.18	2.68	0.412	0.189	-0.16	-0.17	-0.16	-0.101	-0.024	0.02
		qPro-B2-2	Sat_009Satt474	qPro-J-7	Sat_255-Sat_394	MCIM	-0.182	0.38	1.72	-0.277	-0.127	0.132	0.136	0.005	0.016	-0.03	0.144
OC	RIL3613	qOil-D1a-3	Sat_346Satt198	qOil-M-1	Satt567-Satt346	ICIM	-0.025	5.5	1.4	0.197	-0.096	-0.147	-0.383	0.108	0.097	0.223	0.002
		qOil-D1a-1	Sat_332Sat_413	qOil-C1-3	Sat_140-Satt396	MCIM	0.067	0.23	0.46	0.108	0.063	0.009	0.039	-0.056	-0.078	0.013	-0.100
		qOil-M-2	Satt626-Satt536	qOil-F-2	Sat_039-Satt425	MCIM	-0.099	0.34	0.62	0.065	0.01	-0.054	-0.253	0.057	0.048	0.044	0.084
	RIL6013	qOil-H-1	Satt181-Satt434	qOil-D2-2	Sat_194-Sat_001	MCIM	0.073	0.15	0.49	0.007	0.046	0.037	0.072	-0.023	-0.063	-0.064	-0.013
		qOil-C1-2	Sat_367Sat_140	qOil-C2-1	Satt640-Sat_336	MCIM	-0.026	0.04	0.52	0.013	0.147	-0.086	-0.209	0.073	0.015	0.016	0.032

A: The estimated additive by additive epistatic effect

B: Phenotypic variation explained by epistatic QTL

C: Phenotypic variation explained by epistasis × environment (AAE) interactions

D: Epistatic effects by environment interaction. E1: Keshan in 2013; E2: Harbin in 2014; E3: Harbin in 2015; E4: Keshan in 2015; E5: Acheng in 2016; E6: Shuangcheng in 2016; E7: Harbin in 2016; E8: Shuangcheng in 2017

remaining 9 pairwise interaction sites showed positive epistatic effects. For 4 pairwise interaction sites (qPro-M-2~qPro-F-2, qPro-D1b-3~qPro-F-5, qPro-H-3~qPro-G-2 and qPro-D1a-2~qPro-M-3), the PVE for the epistatic QTLs was greater than the PVE for the AAE interaction, indicating that it was strongly impacted by the epistatic effects, whereas for the other nine pairwise interaction sites, the PVE of the epistatic QTLs was lower than that for the AAE interaction, indicating that it is greatly impacted by the environment.

Likewise, we detected five sites of pairwise interaction related to OC by AA and AAE analysis in the two RIL populations (Table 6 and Fig. 3). The AA values ranged from 0.025 to 0.099%, the PVEs for AA ranged from 0.04 to 5.5%, explaining 6.26% of the total variation in OC, and the PVEs for AAE ranged from 0.46 to 1.40%, explaining 3.49% of the total variation in OC. Three pairwise interaction sites, qOil-D1a-3~qOil-M-1, qOil-M-2~qOil-F-2 and qOil-C1-2~qOil-C2-1, showed negative epistatic effects, while the other two (qOil-D1a-1~qOil-C1-3 and qOil-H-1~qOil-D2-2) showed positive epistatic effects. For qOil-D1a-3~qOil-M-1, the PVE for the epistatic effect was greater than that for AAE, indicating that it is greatly impacted by the parents, whereas the reverse was true for the other four pairs epistasis effects QTLs, indicating that it is greatly impacted by the environment.

Fig. 3: Epistatic QTL for protein (blue lines) and oil (red lines) contents in RIL3613 (a) and RIL6013 (b) populations

Overall, for PC and OC combined, among the significantly epistatic QTL pairs that we found, two epistatic effects were due to the interactions of two significant QTLs, seven to the interactions of one significant and one non-significant QTL and the remaining four to the interactions of two non-significant QTLs (Table 4, 5 and 6).

Discussion

RIL populations are homozygous populations in which progeny reliably inherit their parents' traits, generally created by plant breeders as a means to develop new varieties, or to perform QTL mapping (Luo et al. 2015; Warrington et al. 2015). However, the number of polymorphic markers between the parents may be limited, resulting in a low marker density in molecular genetic maps constructed from RILs (Zhang and Wang 2015). To overcome this limitation, plant breeders use multiple-population improvement, a strategy that has been useful in, for instance, rice (Zeng et al. 2017), Arabidopsis thaliana (Bloomer et al. 2014), maize (Li et al. 2014; Pan et al. 2017), soybean (Mao et al. 2013; Kamfwa et al. 2017). However, separate populations may not contain the same QTL markers, making it difficult to accurately estimate the number of common QTLs across multiple genetic backgrounds.

In this study, we used two RIL populations with a common female parent (Dongnong L13) and were able to detect 32 QTLs with overlapping locations in both populations (Fig. 2, 4). The qPro-A1-1 and qOil-A1-1 regions overlapped the qPro-A1-2 region; the qPro-A1-1 region contains a QTL previously found by Mao et al. (2013), while qOil-A1-1 was found to be a hotspot region by Rossi et al. (2013), Brummer et al. (1997), Qi et al. (2011) and Han et al. (2015). In the B1 linkage group, the qPro-B1-1 region contains qPro-B1-2 and qOil-B1-1; the former overlaps with a QTL previously identified by Gai et al. (2007), and also with Seed protein 25-1 (Gai et al. 2007), and the qOil-B1-2 interval contains Seed oil 39-2, identified by Wang et al. (2014b). In the C1 linkage group, the qPro-C1-1 region overlapped qPro-C1-3, identified as a hotspot region found by several previous studies (Orf et al. 1999; Stombaugh et al. 2004; Mao et al. 2013; Wang et al. 2014b). Similarly, the qPro-C2-4 region contains the qPro-C2-5 region, and both are consistent with QTLs identified by Pathan et al. (2013); moreover, numerous QTLs related to seed PC in soybean have been located in the qPro-C2-4 hotspot region (Csanádi et al. 2001; Liang et al. 2010; Pathan et al. 2013; Rossi et al. 2013). In the D1a linkage group, the qPro-D1a-1 and qOil-D1a-2 (5.93–104.27 cM; Sat_413-Sat_160) regions overlapped the qPro-D1a-2 and qPro-D1a-3 regions, the genome is widely located of qOil-D1a-2 and qPro-D1a-1. Several QTLs relevant to soybean protein and oil contents have previously been located in these hotspot regions (Brummer et al. 1997; Csanádi et al. 2001; Specht et al. 2001; Qi et al. 2011; Mao et al. 2013; Wang et al. 2014b; Qi et al. 2014; Han et al. 2015). In addition, qPro-D1a-3 contains Seed protein 40-4 located by Qi et al. (2014). In the D1b linkage group, the qPro-D1b-5 region overlapped the qOil-D1b-1 region detected by Mao et al. (2013) and Qi et al. (2014) and qPro-D1b-3, which controlled PC in both populations, was also found by Qi et al. (2014) in the Charleston and Dongnong 594 soybean strains and can be expressed stably in multiple genetic backgrounds simultaneously. In the F linkage group, the qPro-F-1 region overlapped the qOil-F-1 region and it includes Seed oil 24-4, located by Qi et al. (2011), and is accordant with the QTL identified by Mao et al. (2013). In the G linkage group, the qPro-G-3 region included qOil-G-2

Fig. 4: Genomic region of QTL associated with protein and oil contents in present and previous researches

QTLs shown in red colour and blue colour were identified in RIL3613 and RIL6013 population in this study, respectively; QTLs shown in black were identified in previous studies

and overlapped the Seed protein 20-1 region. In addition, several QTLs related to soybean PC have previously been identified in widely distributed locations of qPro-I-2 (Lu et al. 2013; Rossi et al. 2013; Hacisalihoglu et al. 2018). The qOil-J-1 and qPro-J-1 regions overlapped qPro-J-2, qPro-J-3 and qPro-J-4. Among these, qOil-J-1 is consistent with the results of Mao et al. (2013) and Eskandari et al. (2013). Finally, the qPro-M-4 region overlapped the qPro-M-5 region in the M linkage group.

QTLs can exist in the same chromosome region in different populations simultaneously, which can to some extent allow the improvement of multiple traits at the same time. Here, we compared newly identified QTLs with those known from previous studies of strains with different genetic backgrounds to improve the accuracy and versatility of these QTLs.

Some specific QTLs were identified in only one of our two populations. In this study, 41 QTLs located on 11 chromosomes (A1, A2, B2, C1, C2, D2, F, G, H, I and O) were found to have no overlapping region in the two mapping groups (Fig. 2 and 4). Most of the QTLs we found, with the exception of qPro-A1-3, qPro-A2-2, qOil-C1-1, qOil-C2-3 and qOil-H-1, were already known from previous studies. Only some are stable in different genetic backgrounds—such as qPro-I-1, qPro-C2-3, qOil-B2-1, qOil-C2-2, qOil-D2-1 and qOil-I-1, which are in known hotspot regions—which underlines the potential importance of the influence of specific QTLs in breeding.

Beside confirming various QTLs already found to be as associated with soybean protein and oil contents in previous researches (as discussed above), we also identified 10 previously unknown QTLs in the RIL3613 population and seven QTLs in the RIL6013 population that are associated with one or both of these traits.

Neglecting the presence of epistasis impairs the ability to recognize QTLs and reduces the efficiency of MAS (Palomeque et al. 2010; Korir et al. 2011; Qi et al. 2017). We therefore mapped the epistatic effects (AA) and epistasis by environment interaction effect (AAE) for PC and OC using ICIM and MCIM models for two RIL populations in eignt environments. Overall, we detected 13 and five epistatic QTL pairs for PC and OC, respectively, in linkage groups B2, C1, C2, D1a, D1b, D2, F, G, H, J, L, M and N. Traits are affected not only by main effect QTLs but also by the interactions among loci (Ding et al. 2014; Jannink 2007; Tan et al. 2018); thus, epistatic effects are a significant factor for complex traits, such as PC and OC. In the present study, the multi-environment joint analysis method identified two pairs of epistatic QTLs that occur between significantly additive QTLs, as well as 6 significant additive effects QTLs, that participate in epistatic and environmental interactions, interact with other QTLs, and increase the phenotypic variation of the epistasis effect, the overall phenotypic variation and the MAS efficiency, as indicated by the phenotypic variation explained (PVE) value of significant additive effect (Fig. 3, Table 4, 5 and 6 underline). The other five pairs of epistatic QTLs are linked by non-significant additive QTLs, which indicates that QTL can not only directly affect phenotypic expression, but also affect the expressed traits through interactions with other loci; this knowledge can be used to improve the efficacy of QTL detection, which is related to the general genetic status of quantitative traits (Li et al. 2014; Teng et al. 2017). Four of these pairs of epistatic QTLs, qPro-D1a-2~qPro-M-3, qOil-D1a-3~qOil-M-1, qPro-D1a-4~qPro-M-4 and qPro-D1a-4~qPro-M-1, involve QTLs located in linkage groups D1a and M; two other pairs of epistatic QTLs, qPro-M-2~qPro-F-2 and qOil-M-2~qOil-F-2, are between QTLs in linkage groups M and F; and qPro-J-6~qPro-J-4 and qPro-J-6~qPro-J-5, are between QTLs in the same linkage group. qPro-D1a-4 and qPro-J-6 are stable loci whose epistatic interaction has been repeatedly identified, and it seems plausible that they may contain genes regulating PC in soybean seeds. The above six pairs of epistatic QTL regions all overlap to some extent, indicating that a QTL controlling one trait may produce multiple epistatic effects in different environments.

Many studies have shown that the PC and OC of soybean seeds can be affected by common markers, but there have been relatively few studies showing the influence from overlap between common epistatic interaction regions (Brummer et al. 1997; Csanádi et al. 2001; Lee et al. 2019). Here, we found that the overlapping qPro-D1a-2~qPro-M-3 (marker interval Sat_346-Satt515~Satt567-Satt697) and qOil-D1a-3~qOil-M-1 (marker interval Sat_346-Satt198~Satt567-Satt346) regions jointly control soybean seed PC and OC, as do qPro-M-2~qPro-F-2 (marker interval Sat_389-Satt697~Satt030-Sat_262) and qOil-M-2~qOil-F-2 (marker interval Satt626-Satt536~Sat_039-Satt425) (Fig. 3 and Table 6). These results indicate that epistatic interaction plays a major role in the accumulation of PC and OC in soybean seed and must be taken into consideration in investigating the genetic bases of these two traits.

Epistatic effects and environmental factors play major roles to formation in complex traits (Allard 1996; Karikari et al. 2019). Soybean seed protein and oil content QTLs have genetic specificity and environmental sensitivity (Wang et al. 2015b) and can thus be identified by analysis of AA and AAE QTLs. A low PVE for AAE indicates an epistatic effect is non-essentially affected by the environment, and thus a QTL that can be stably expressed, whereas a high PVE for AAE indicates a highly environmentally sensitive QTL. In this study, the PVEs for 5 epistatic QTL pairs, qPro-M-2~qPro-F-2, qPro-D1b-3~qPro-F-5, qPro-H-3~qPro-G-2, qPro-D1a-2~qPro-M-3 and qOil-D1a-3~qOil-M-1, were greater than the PVEs for environmental interaction, indicative of stable inheritance in different environments, whereas the remaining epistatic QTL pairs are environmentally sensitive and only expressed in particular environments (Fig. 3 and Table 6). In MAS breeding strategies for seed protein and oil traits, it is important not to merely consider the additive and epistatic effect QTLs, and additive × environment (AE) and epistasis × environment (AAE) interaction effect QTLs must also be considered for a specific environment. Stabilizing effect QTLs with weak or no interaction with the environment, stable genetic bases and high degrees of variation should be selected.

Conclusion

We detected 50 PC and 23 OC additive-effect QTLs and 13 PC and 5 OC epistatic-effect QTL pairs in two soybean populations. Of these, 12 QTLs were in previously known hotspot regions and 17 QTLs were newly identified, giving these results theoretical and practical significance for future MAS initiatives.

Acknowledgements

The authors gratefully acknowledge the financial support for this study provided by the grants from the National Key Research and Development Program of China (2017YFD0101303) to W-X.L., Project of Research and Development on Applied Technology of Harbin in Heilongjiang Province [China] (2017RAXXJ019) to H.N.

References

Allard RW (1996). Genetic basis of the evolution of adaptedness in plants. Euphytica 92:1‒11

Bloomer RH, AM Lloyd, V Symonds (2014). The genetic architecture of constitutive and induced trichome density in two new recombinant inbred line populations of Arabidopsis thaliana: Phenotypic plasticity, epistasis, and bidirectional leaf damage response. BMC Plant Biol 14; Article 119

Brummer EC, GL Graef, J Orf, JR Wilcox, RC Shoemaker (1997). Mapping QTL for seed protein and oil content in eight soybean populations. Crop Sci 37:370‒378

Chiari L, ND Piovesan, LK Naoe, IC José, JMS Viana, MA Moreira, EG Barros (2004). Genetic parameters relating isoflavone and protein content in soybean seeds. Euphytica 138:55‒60

Csanádi G, J Vollmann, G Stift, T Lelley (2001). Seed quality QTLs identified in a molecular map of early maturing soybean. Theor Appl Genet 103:912‒919

Ding G, Z Zhao, L Wang, D Zhang, L Shi, F Xu (2014). Identification and multiple comparisons of QTL and epistatic interaction conferring high yield under boron and phosphorus deprivation in Brassica napus. Euphytica 198:337‒351

Eskandari M, ER Cober, I Rajcan (2013). Genetic control of soybean seed oil: II. QTL and genes that increase oil concentration without decreasing protein or with increased seed yield. Theor Appl Genet 126:1677‒1687

Gai J, Y Wang, X Wu, SA Chen (2007). Comparative study on segregation analysis and QTL mapping of quantitative traits in plants with a case in soybean. Front Agric Chin 1:1‒7

Hacisalihoglu G, AL Burton, JL Gustin, S Eker, S Asikli, EH Heybet, L Ozturk, I Cakmak, A Yazici, KO Burkey, J Orf, AM Settles (2018). Quantitative trait loci associated with soybean seed weight and composition under different phosphorus levels. J Integr Plant Biol 60:232‒241

Han Y, W Teng, Y Wang, X Zhao, L Wu, D Li (2015). Unconditional and conditional QTL underlying the genetic interrelationships between soybean seed isoflavone, and protein or oil contents. Plant Breed 134:300‒309

Hou M, Z Qi, X Han, D Xin, H Jiang, C Liu, Q Wu, L Sui, G Hu, Q Chen (2014). QTL mapping and interaction analysis of seed protein content and oil content in soybean. Sci Agric Sin 47:2680‒2689

Jannink J (2007). Identifying quantitative trait locus by genetic background interactions in association studies. Genetics 176:553‒561

Kamfwa K, D Zhao, JD Kelly, KA Cichy (2017). Transcriptome analysis of two recombinant inbred lines of common bean contrasting for symbiotic nitrogen fixation. PLoS One 12; Article e0172141

Karikari B, S Li, J Bhat, Y Cao, J Kong, J Yang, T Zhao (2019). Genome-wide detection of major and epistatic effect QTLs for seed orotein and oil content in soybean under multiple environments using high-density bin map. Intl J Mol Sci 20:979‒999

Korir PC, B Qi, Y Wang, T Zhao, D Yu, S Chen, J Gai (2011). A study on relative importance of additive, epistasis and unmapped QTL for Aluminium tolerance at seedling stage in soybean. Plant Breed 130:551‒562

Lee S, K Van, M Sung, R Nelson, MAR Mian (2019). Genome-wide association study of seed protein, oil and amino acid contents in soybean from maturity groups i to iv. Theor Appl Genet 132:1639‒1659

Li H, G Ye, J Wang (2006). A Modified Algorithm for the Improvement of Composite Interval Mapping. Genetics 175:361‒374

Li K, J Yan, J Li, X Yang (2014). Genetic architecture of rind penetrometer resistance in two maize recombinant inbred line populations. BMC Plant Biol 14:1471‒1482

Liang HZ, YL Yu, SF Wang, Y Liang, TF Wang, YL Wei, PT Gong, XY Liu, XJ Fang, MC Zhang (2010). QTL mapping of isoflavone, oil and protein contents in soybean (Glycine max L. Merr.). Agric Sci Chin 9:1108‒1116

Lu W, Z Wen, H Li, D Yuan, J Li, H Zhang, Z Huang, S Cui, W Du (2013). Identification of the quantitative trait loci (QTL) underlying water soluble protein content in soybean. Theor Appl Genet 126:425‒433

Luo J, SA Jobling, A Millar, MK Morell, Z Li (2015). Allelic effects on starch structure and properties of six starch biosynthetic genes in a rice recombinant inbred line population. Rice 8:15-28

Mao T, Z Jiang, Y Han, W Teng, X Zhao, W Li (2013). Identification of quantitative trait loci underlying seed protein and oil contents of soybean across multi-genetic backgrounds and environments. Plant Breed 132:630‒641

McCouch SR, YG Cho, M Yano, E Paul, M Blinstrub (1997). Report on QTL nomenclature. Rice Genet Newsl 14:11‒13

Ning H, J Yuan, Q Dong, W Li, H Xue, Y Wang, Y Tian, WX Li (2018). Identification of QTLs related to the vertical distribution and seed-set of pod number in soybean [Glycine max (L.) Merri]. PLoS One 13; Article e0195830

Orf JH, K Chase, T Jarvik, LM Mansur, KG Lark (1999). Genetics of soybean agronomic traits: I. Comparison of three related recombinant inbred populations. Crop Sci 39:1642‒1651

Palomeque L, L Liu, W Li, BR Hedges, ER Cober, MP Smid, L Lukens, I Rajcan (2010). Validation of mega-environment universal and specific QTL associated with seed yield and agronomic traits in soybeans. Theor Appl Genet 120:997‒1003

Pathan SM, T Vuong, K Clark, JD Lee, DA Sleper (2013). Genetic Mapping and Confirmation of Quantitative Trait Loci for Seed Protein and Oil Contents and Seed Weight in Soybean. Crop Sci 53:765‒774

Pan Q, Y Xu, K Li, Y Peng, W Zhan, W Li, L Li, J Yan (2017). The genetic basis of plant architecture in 10 maize recombinant inbred line populations. Plant Physiol 175:858-873

Qi Z, X Zhang, H Qi, D Xin, X Han, H Jiang, Z Zhang, J Zhang, R Zhu, Z Hu, C Liu, X Wu, Q Chen, Z Yin, C Daidi (2017). Identification and validation of major QTLs and epistatic interactions for seed oil content in soybeans under multiple environments based on a high-density map. Euphytica 213:162‒175

Qi Z, M Hou, X Han, C Liu, H Jiang, D Xin, G Hu, Q Chen (2014). Identification of quantitative trait loci (QTLs) for seed protein concentration in soybean and analysis for additive effects and epistatic effects of QTLs under multiple environments. Plant Breed 133:499‒507

Qi Z, Q Wu, X Han, Y Sun, X Du, C Liu, H Jiang, G Hu, Q Chen (2011). Soybean oil content QTL mapping and integrating with meta-analysis method for mining genes. Euphytica 179:499‒514

Rossi ME, JH Orf, L Liu, Z Dong, I Rajcan (2013). Genetic basis of soybean adaptation to North American vs. Asian mega-environments in two independent populations from Canadian × Chinese crosses. Theor Appl Genet 126:1809‒1823

Specht JE, K Chase, M Macrander, GL Graef, J Chung, JP Markwell, JHO Germann, KG Lark (2001). Soybean response to water: A QTL analysis of drought tolerance. Crop Sci 41:493‒509

Stombaugh SK, JH Orf, HG Jung, K Chase, KG Lark, DA Somers (2004). Quantitative Trait Loci Associated with Cell Wall Polysaccharides in Soybean Seed. Crop Sci 44:2101‒2106

Tan R, B Serven, PJ Collins, Z Zhang, Z Wen, JF Boyse, C Gu, MI Chilvers, BW Diers, D Wang (2018). QTL mapping and epistatic interaction analysis of field resistance to sudden death syndrome (Fusarium virguliforme) in soybean. Theor Appl Genet 131:1729‒1740

Teng W, B Zhang, Q Zhang (2017). Identification of quantitative trait loci underlying seed oil content of soybean including main, epistatic and QTL × environment effects in different regions of Northeast China. Crop Past Sci 68:625‒631

Wang X, G Jiang, M Green, RA Scott, Q Song, DL Hyten, PB Cregan (2014a). Identification and validation of quantitative trait loci for seed yield, oil and protein contents in two recombinant inbred line populations of soybean. Mol Genet Genomics 289:935‒949

Wang X, G Jiang, M Green, RA Scott, DL Hyten, PB Cregan (2014b). Quantitative trait locus analysis of unsaturated fatty acids in a recombinant inbred population of soybean. Mol Breed 33:281‒296

Wang Y, Y Han, X Zhao, Y Li, W Teng, D Li, Y Zhan, W Li (2015a). Mapping isoflavone QTL with main, epistatic and QTL × environment effects in recombinant inbred lines of soybean. PLoS One 10; Article e0118447

Wang J, P Chen, D Wang, G Shannon, A Zeng, M Orazaly, C Wu (2015b). Identification and mapping of stable QTL for protein content in soybean seeds. Mol Breed 35:92‒101

Warrington CV, H Abdel-Haleem, DL Hyten, PB Cregan, JH Orf, AS Killam, N Bajjalieh, Z Li, HR Boerma (2015). QTL for seed protein and amino acids in the Benning × Danbaekkong soybean population. Theor Appl Genet 128:839‒850

Yang J, J Zhu (2005). Methods for predicting superior genotypes under multiple environments based on QTL effects. Theor Appl Genet 110:1268‒1274

Yang J, C Hu, C Hu, R Yu, Z Xia, X Ye, J Zhu (2008). QTL Network: Mapping and visualizing genetic architecture of complex traits in experimental populations. Bioinformatics 24:721‒723

Zeng Y, J Shi, Z Ji, Z Wen, Y Liang, C Yang (2017). Combination of twelve alleles at six quantitative trait loci determines grain weight in rice. PLoS One 12; Article e0181588

Zhang H, H Wang (2015). QTL mapping for traits related to P-deficient tolerance using three related RIL populations in wheat. Euphytica 203:505‒520